function spatialAnalysis(cellType)
%  FUNCTION SPATIALANALYSIS(CELLTYPE) looks at the relationship between the
%  connection probability and the location of the cells.

% load R, Theta, and connectionBoolean
load(['C:\Data\Taro\ANALYSIS\multi_rec_project\connectivity\cellDistributionAnalysis\',cellType,'Scatter.mat']);

connectionIndex = find(connectionBoolean == 1);
noconnectionIndex = find(connectionBoolean == 0);

ybasedDist(Theta, R, connectionIndex, noconnectionIndex)

% RThetabasedDist(R, Theta, connectionIndex, noconnectionIndex)

distEdge = [0 150 300];
angleEdge = [-pi -pi/2 pi/2, pi];
[Con180 Uncon180] = R4divideDist(R, Theta, connectionIndex, noconnectionIndex, angleEdge, 1, distEdge);
hpatch = findobj(gcf,'Type','patch');
set(hpatch,'UserData',{Con180,Uncon180});

angleEdge = [-pi -pi/2 0 pi/2 pi];
[Con4pi Uncon4pi] = R4divideDist(R, Theta, connectionIndex, noconnectionIndex, angleEdge, 0, distEdge);
hpatch = findobj(gcf,'Type','patch');
set(hpatch,'UserData',{Con4pi,Uncon4pi});

angleEdge = [-pi -pi*3/4 -pi/4 pi/4 pi*3/4 pi];
[ConTilt UnconTilt] = R4divideDist(R, Theta, connectionIndex, noconnectionIndex, angleEdge, 1, distEdge);
hpatch = findobj(gcf,'Type','patch');
set(hpatch,'UserData',{ConTilt,UnconTilt});

% save(['C:\Data\Taro\ANALYSIS\multi_rec_project\connectivity\cellDistributionAnalysis\bullseyeSummaryData',cellType,'.mat'],...
%     'Con180','Uncon180','Con4pi','Uncon4pi','ConTilt','UnconTilt');
end
function [xVal, yVal] =  ConTable2XY(ContingencyTable)

xVal = [];
yVal = [];
for kk = 1:size(ContingencyTable, 2)
    xVal = [xVal,ones(1,ContingencyTable(1,kk) + ContingencyTable(2,kk)) * kk];
    yVal = [yVal,zeros(1,ContingencyTable(2,kk)),ones(1,ContingencyTable(1,kk))];
end
end

function ybasedDist(Theta, R, connectionIndex, noconnectionIndex)
% first, devide the area based on y position.
[x, y] = pol2cart(Theta + pi/2, R); % R is defined as the angle from the radial axis. Rotate pi/2 so y axis represents the radial axis.

yCon = y(connectionIndex); % y coordinates of connected postsynaptic cells.
yUncon = y(noconnectionIndex); % y coordinates of unconnected postsynaptic cells.

% classify yCon and yUncon based on their y position.
edge = [-inf -100 -50 0 50 100 inf];
[nCon, binCon] = histc(yCon, edge);
[nUncon, binCon] = histc(yUncon, edge);

% plot the distribution of connectivity.
connectRate = nCon ./ (nCon + nUncon);
figure;
barh(connectRate(1:length(connectRate) - 1));
for k = 1:length(edge)-1
    xl{k} = [num2str(edge(k)), ' to ', num2str(edge(k + 1))];
end
set(gca,'YTickLabel', xl)
% saveas(gcf, ['C:\Data\Taro\ANALYSIS\multi_rec_project\connectivity\cellDistributionAnalysis\',cellType,'ydist.fig'])
xyContingencyTable = [nCon;nUncon];

% independence test. H0: connectivity does not depend on the y position.
[xVal, yVal] = ConTable2XY([nCon;nUncon]);
[sig,p,con] = FisherExactTest(xVal, yVal)
disp('warning: sig =0 means dependency...')
end

function RThetabasedDist(R, Theta, connectionIndex, noconnectionIndex)
% second, divide the area based on r and theta.
rCon = R(connectionIndex); % r vals of connected postsynaptic cells.
rUncon = R(noconnectionIndex); % r vals of unconnected postsynaptic cells.
thetaCon = Theta(connectionIndex);
thetaUncon = Theta(noconnectionIndex);

% classify the cells based on, 1. connected vs unconnected, 2. angle, 3.
% radius.
thetaConUpperIndex = intersect(find(thetaCon > -pi / 2), find(thetaCon < pi / 2));
nUpCon = histc(rCon(thetaConUpperIndex),[0 50 100 inf]);

thetaConLowerIndex = union(find(thetaCon < -pi / 2), find(thetaCon > pi / 2));
nLowCon = histc(rCon(thetaConLowerIndex),[0 50 100 inf]);

thetaUnconUpperIndex = intersect(find(thetaUncon > -pi / 2), find(thetaUncon < pi / 2));
nUpUncon = histc(rUncon(thetaUnconUpperIndex),[0 50 100 inf]);

thetaUnconLowerIndex = union(find(thetaUncon < -pi / 2), find(thetaUncon > pi / 2));
nLowUncon = histc(rUncon(thetaUnconLowerIndex),[0 50 100 inf]);

% Fisher exact test for contingency.
ContingencyTable = [fliplr(nUpCon), nLowCon; fliplr(nUpUncon), nLowUncon];
ContingencyTable = ContingencyTable(:,2:size(ContingencyTable,2)-1)
[xVal, yVal] = ConTable2XY(ContingencyTable);
[sig,p,con] = FisherExactTest(xVal, yVal)
disp('warning: sig =0 means dependency...')

figure;
bullsConnectRate = [nUpCon ./ (nUpCon + nUpUncon); nLowCon ./ (nLowCon + nLowUncon)];
bullsConnectRate = bullsConnectRate(:,1:size(bullsConnectRate,2) - 1);
bullseye(bullsConnectRate ,'n',100,'rho',[0 3],'tht',[90 450])
% saveas(gcf, ['C:\Data\Taro\ANALYSIS\multi_rec_project\connectivity\cellDistributionAnalysis\',cellType,'bullseye.fig'])

% chisquare in ~100 um.
Contingency100 = ContingencyTable(:,size(ContingencyTable, 2)/2-1:size(ContingencyTable, 2)/2+2);
Contingency100 = [sum(Contingency100(:,1:2),2),sum(Contingency100(:,3:4),2)]
chisquarecont(Contingency100)
end

function [Con Uncon] = R4divideDist(R, Theta, connectionIndex, noconnectionIndex, angleEdge,fuseFlag,distEdge)
% second, divide the area based on r and theta.
rCon = R(connectionIndex); % r vals of connected postsynaptic cells.
rUncon = R(noconnectionIndex); % r vals of unconnected postsynaptic cells.
thetaCon = Theta(connectionIndex);
thetaUncon = Theta(noconnectionIndex);

% classify the cells based on, 1. connected vs unconnected, 2. angle, 3.
% radius.



for k = 1:length(angleEdge) - 1
    ThetaConnectionIndex{k} = intersect(find(thetaCon > angleEdge(k)), find(thetaCon <= angleEdge(k + 1)));
    Con{k,1} = histc(rCon(ThetaConnectionIndex{k}),distEdge);

    ThetaUnconnectionIndex{k} = intersect(find(thetaUncon > angleEdge(k)), find(thetaUncon <= angleEdge(k + 1)));
    Uncon{k,1} = histc(rUncon(ThetaUnconnectionIndex{k}),distEdge);
end

Con = cell2mat(Con);
Uncon = cell2mat(Uncon);
if fuseFlag == 1
    areaNum = size(Con, 1);
    Con(1,:) = Con(1,:) + Con(areaNum,:);
    Con(areaNum,:) = [];
    Uncon(1,:) = Uncon(1,:) + Uncon(areaNum,:);
    Uncon(areaNum,:) = [];
    tht = [-180/(areaNum-1) 360 - 180 / (areaNum-1)];
else
    tht = [0 360];
end
Con
Uncon

figure;



bullsConnectRate = Con ./ (Uncon + Con)
bullsConnectRate = bullsConnectRate(:,1:size(bullsConnectRate,2) - 1);
bullseye(bullsConnectRate ,'n',100,'rho',[0 3],'tht',tht)

% % Fisher exact test for contingency.
% ContingencyTable = [fliplr(nUpCon), nLowCon; fliplr(nUpUncon), nLowUncon];
% ContingencyTable = ContingencyTable(:,2:size(ContingencyTable,2)-1)
% [xVal, yVal] = ConTable2XY(ContingencyTable);
% [sig,p,con] = FisherExactTest(xVal, yVal)
% disp('warning: sig =0 means dependency...')
% 
% % chisquare in ~100 um.
% Contingency100 = ContingencyTable(:,size(ContingencyTable, 2)/2-1:size(ContingencyTable, 2)/2+2);
% Contingency100 = [sum(Contingency100(:,1:2),2),sum(Contingency100(:,3:4),2)]
% chisquarecont(Contingency100)
end
